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SUMMARY 


As  part  of  the  study  of  the  nonlinear  response  of  structural 
components  we  investigate  in  this  report  a  single  degree-of-freedom  non¬ 
linear  response  equation  usually  known  as  Duffing’s  equation.  Previously  the 
response  curves  have  been  investigated  analytically  by  many  authors  but 
there  seems  to  have  been  very  little  numerical  work  done  in  the  nonlinear 
cases.  It  is  shown  that  a  numerical  solution  is  feasible  and  that  some  of  the 
standard  analytical  approximations  are  not  accurate  in  the  peak  region  of  the 
response  curve. 


RESUME 


Dans  le  cadre  d’une  etude  de  la  reponse  non-lineaire  d’elements 
structuraux,  les  auteurs  analysent  une  equation  a  un  seul  degre  de  liberte 
relative  a  la  reponse  non-lineaire,  connue  generalement  sous  le  nom  d’equa- 
tion  de  Duffing.  Les  courbes  de  reponse  ont  deja  ete  analysees  par  de  nom- 
breux  auteurs,  mais  il  semble  qu’on  ait  fait  tres  peu  de  travaux  de  nature 
numerique  sur  les  cas  non-lineaires.  Les  auteurs  demontrent  qu’une  solution 
numerique  est  possible  et  que  certaines  des  approximations  analytiques 
standard  ne  sont  pas  exactes  en  ce  qui  conceme  le  pic  de  la  courbe  de 
reponse. 
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1.0  INTRODUCTION 
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In  recent  years  there  have  been  rapid  developments  in 
structural  analysis  techniques  and  related  computer  codes.  The 
analysis  of  linear  response  of  structural  components  has  been  well 
studied  and  numerous  numerical  time  marching  techniques  have  been 
developed  for  finite  element  structural  analysis  programs.  Some  of  the 
more  commonly  used  time  integration  techniques  are  the  explicit  central 
difference  scheme  and  Houbolt's,  Wilson's  and  Newmark’s  methods 
(Ref.  1).  In  order  for  any  algorithm  to  be  useful  in  practice,  it  must 
not  lead  to  a  divergent  solution.  Except  for  the  explicit  scheme,  the 
other  three  methods  are  unconditionally  stable.  A  study  of  the 
accuracies  shows  that  all  three  schemes  suffer  from  a  period  elongation 
which  increases  with  increase  in  time  step.  In  addition,  Houbolt's  and 
Wilson's  methods  have  artifical  damping  resulting  in  the  higher  modes 
being  rapidly  damped  out  in  a  system  with  a  large  number  of  degrees  of 
freedom. 

In  non-linear  problems  established  methods  for  investigating 
stability  and  decay  as  in  linear  responses  are  not  possible.  Various 
formulations  have  been  proposed  and  used  with  varying  degrees  of 
success.  The  four  methods  described  in  [1]  for  linear  problems  can 
readily  be  used  for  the  nonlinear  case.  The  choice  of  any  particular 
algorithm  depends  to  a  large  extent  on  the  type  of  problem  to  be 
solved. 

The  use  of  time  marching  techniques  in  flutter  studies  was 
first  reported  by  Ballhaus  and  Goorjian  [2]  who  carried  out  the 
aeroelastic  response  of  a  NACA  64A006  airfoil  with  a  single 


degree-of-f reedom  in  pitch  at  transonic  speeds.  Extensions  of  this 
procedure  to  two-  and  three-degree-of-f reedom  have  been  reported  by 
Rizzetta  [3]  and  Yang  and  Chen  [4],  Only  a  linear  response  was  treated 
by  these  authors.  In  flutter  studies,  there  are  potentially  many 
sources  of  non-  linearities  present,  some  of  which  may  give  rise  to 
"limit-amplitude''  oscillations.  The  nonlinearities  may  arise  in  many 
possible  ways,  but  the  most  commonly  encountered  ones  are  those  having 
structural  or  aerodynamic  origin. 

In  aeroelastic  applications  unconditionally  stable  methods 
like  Houbolt's  for  linear  problems  is  not  very  critical,  since  the 
number  of  modes  considered  is  usually  not  very  large.  Conditionally 
stable  schemes  can  then  be  used,  provided  a  time  step  is  chosen  to 
ensure  the  highest  mode  does  not  diverge.  Also,  it  is  important  that 
amplitude  decay  and  period  elongation  should  be  as  small  as  possible. 

In  the  present  study  higher  order  implicit  finite  difference  schemes 
with  greater  accuracies  than  Houbolt's  method  are  discussed.  To  study 
the  stability  and  accuracies  of  the  higher  order  schemes,  a  linear 
system  is  considered. 

In  the  application  to  the  nonlinear  case  the  restoring  force 
is  assumed  to  be  given  by  that  of  a  cubic  spring.  The  equation  to  be 
solved  is  essentially  Duffing's  equation  [5],  In  aeroelastic  problems 
a  cubic  restoring  force  corresponding  to  a  hardening  spring  arises  when 
a  thin  wing  or  propeller  is  subjected  to  increasingly  large  amplitudes 
in  torsion. 

The  numerical  solution  to  the  nonlinear  Duffing's  equation 
offered  some  interesting  studies.  For  instance,  it  was  not  known  how 
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a  multivalued  solution  could  be  computed  and  whether  the  jump 
phenomenon  from  one  branch  of  the  response  curve  to  another  could  be 
forecast  numerically.  What  would  affect  such  a  jump  and  would  the  jump 
be  consistent?  Also  would  the  unstable  part  of  the  response  curve 
cause  numerical  difficulties?  Are  the  approximate  analytic  response 
curves  accurate? 

We  present  here  solutions  which  answer  the  above  questions 
and  hopefully  this  will  clear  the  way  to  study  higher  order  degrees-of- 
freedom  systems  and  to  more  practical  applications. 


2.0  THE  RESPONSE  EQUATIONS 
2.1  The  Linear  Equation 

In  order  to  study  our  numerical  schemes  we  first  consider  the 
linear  equation  for  damped  oscillations  i.e. 

x  +  cx  +  kx  =  F  sin  wt 


Its  solution  is 

— r*  t*  /  9  O  n 

x  =  e  [A  sin  /(I -£  )  a)  t  +  A, cos  /(I-5  )  gj  t] 

2  n  J  n 

+  A^  sin(<ot-(J> ) 
where 


a3  =  +  A4  sin(j> 


A2  =  Lxo  +  A3  “  A^ojcos^]  // (l-£  ) 

a,  =  [U-n2)2  +  (2?n)2]  * 


Z  -  ;  u  =  /k  ;  SI  »  — 

2/k  n 

2zn 


i-fi 


(D 

n 


tam}> 
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and  Xq,  Xq  are  the  initial  displacement  and  velocity  respectively.  It 
can  be  seen  that  transients  represented  by  the  first  term  for  x  die 
down  and  leave  the  particular  integral  A^  sin(u)t-if> ).  Thus  the  initial 
conditions  x^,  x^  are  finally  unimportant  in  this  linear  case.  Only  in 
the  nonlinear  system  will  they  be  critical  as  shown  by  the  numerical 
solutions . 

2.2  The  Nonlinear  Equation 

Many  articles  have  been  written  on  the  equation  which 
represents  forced  oscillations  of  systems  with  a  nonlinear  restoring 
force.  The  equation  can  be  written 

x  +  cx  +  f(x)  =  F  sin  u>t 

for  c  >  0.  Most  of  the  studies  have  centred  on  displacements  large 

enough  to  include  linear  and  cubic  terms  in  f(x)  i.e.  the  equation 

"  3 

x  +  cx  +  kx  +  8x  =  F  sin  u>t  (1) 

is  normally  considered  and  this  is  the  form  we  will  concentrate  upon  in 
this  study.  This  equation  is  normally  referred  to  as  Duffing's 
equation,  since  Duffing  [5]  was  the  first  to  make  significant  progress 
in  studying  it. 

Duffing  in  his  study  neglected  the  damping  force  cx  and,  by 

an  expansion  method,  found  that  to  first  order  the  equation 

2  3  2  F 

a>  =  k  +  I  0AZ  -  j  (2) 

related  oj  to  the  amplitude  A  of  the  oscillation.  The  equation  (2)  is 
an  interesting  relation  in  that  for  certain  values  of  the  parameters 
there  are  three  values  of  A  associated  with  one  value  of  w  as 
Illustrated  in  Fig.  1.  This  gave  rise  to  speculation  that  a  jump 
phenomenon  might  exist  in  which  the  solution  would  suddently  jump  from 


one  branch  5-4-6  say  In  Fig.  1  to  the  branch  1-2-3  with  the  branch  3-6 
being  unstable.  A  hysteresis  effect  may  possibly  be  observed  if  the 
amplitude  follows  the  path  5-4-6-2-1  for  increasing  go  and  the  path 
1-2-3-4-5  for  decreasing  go.  Figure  1  shows  the  situation  for  0  >  0 
(hard  spring  i.e.  the  stiffness  increases  with  displacement).  For 
0  <  0  (soft  spring)  the  curve  tilts  over  to  the  left.  Experiments 
confirming  this  hysteresis  effect  are  reported  in  [6], 

A  first  order  expansion  can  also  be  developed  for  equation 


(1)  which  includes  the  damping  term  cx.  This  yields  [see  7] 


2  2  + 

2  i  A  3  a*2  +  F  Ti  c  A%* 
GO  =1+-0A  ±  -  [1 - =-J 

F 


(3) 


4  ~  A 

where  we  have  assumed  that  the  linear  coefficient,  k,  is  1.  From  (3) 
we  observe  that 

I? 

(4) 


Al  <  I 

c 


where  F  is  assumed  positive.  A  different  expansion  using  the  method  of 
Kryloff  and  Bogoliuboff  [8]  yields 


2  2  2  + 
2  ,  .  3  ..2  c  .  F  r,  c  ^  3  2  c  xiT 

(o  =  1+t-BA  ~  T~  ±  T  U"  [1  +  x  BA  --HJ 


(5) 


from  which  an  upper  limit  for  A  can  be  found  namely 

2  /  2  2  i 

_ .  r  i  _  iLJ  +  r  _A_  r  i  _  ^  M  - 

30 


-wd-b)  +-  [Ad-r) 


r] 


(6) 


90  30c 

in  which  the  +  sign  will  be  used  for  0  >  0. 

Comparisons  of  the  maximum  amplitude  given  by  eqs  (4)  and  (6) 


differ  significantly  as  shown  later  where,  in  the  case  studied,  the 
estimate  from  eq  (4)  is  grossly  misleading. 


It  should  be  stressed  that  these  approximations  do  not  yiexd 
any  information  about  the  effect  of  the  initial  conditions  on  the  final 
solution  (after  the  transients  have  died  down).  In  other  words  the 
effect  of  the  initial  conditions  on  which  branch  4-6  or  2-3  the 
solution  lies  is  not  clear  from  the  analytic  studies.  Only  numerically 
can  we  assess  the  effect  of  the  initial  conditions. 


3.0  INTEGRATION  METHODS 

3.1  Central  Difference  Method 

If  we  replace  first  and  second  derivatives  in  (1)  at  time  t 

by  central  difference  expressions  we  obtain 

x  -  2x  +  x  .  x  -  x  i 

n+1  n  n-1  .  n+1  n-1  .  ,  ,  „  3  „  , 

- - -  +  c  -  -  +  kx  +  3x  =  F  sin  wt 

.2  2At  n  n  n 

At 

from  which  x^^  can  be  found.  Thus  the  step  from  time  t  to  t  +  At  can 

4 

be  made  with  accuracy  0(At  )  since  the  central  difference  formulas  are 
2 

of  0(At  )  for  the  derivatives. 

However  this  explicit  method  is  unstable  unless  At  is  less 
than  a  critical  value  and  so  in  practice  the  method  is  not  normally 
used  [1],  Implicit  methods  given  below  have  better  stability  and/or 
accuracy. 

3.2  Houbolt 's  Method 

In  this  case  we  replace  the  derivatives  at  time  t  +  At  with 
backward  difference  formulas  using  values  at  three  previous  points. 


[  2x  -  5x  +  4x  -  x  _]  +  0(At2) 
2  1  n+1  n  n-1  n-2J 


x  =*  — - —  [llx  -  18x  +  9x  -  2x  ]  +  O(At^) 

n  6At  n+1  n  n-1  n-2 


-.1 


These  formulas  and  others  are  found  in  convenient  form  in  Appendix  III 
of  Kopal's  book  {9], 


On  substituting  these  formulas  Into  (1),  written  at  t  +  At, 


we  obtain  an  Implicit  formula  for  xn+^»  namely 


fix  ,  +  Bx  .  «  R 
n+1  n+1 


(7) 


with  accuracy  0(At  )  on  each  step,  where 

B 


__2_  +  JLL£.  +  k 

it2  64t 


F  sin  o>t 


5x  -4x  ,+x  « 

+  _J1 - S  l . -~2  +  -|_  [i8x  -9x  +2x  ] 

6At  L  n  n-1  n-2J 


n+1  .  2 

At 


The  cubic  equation  (7)  Is  solved  by  exact  methods  (see  for  example 
[10])  using  Fortran's  complex  arithmetic  ability.  Since  B  Is  positive 
(c  and  k  being  positive)  it  can  be  seen  from  (7)  that  there  is  only  one 
real  root  if  3  is  positive.  If  3  is  negative  (but  small)  there  will  be 
three  real  roots  since  in  (7)  the  left  hand  side  tends  to  +»  for 


n+1 


and  vice  versa  for  xn+^  +  '+co  whereas  for  typical  values  of 


xn+^  the  predominant  slope  of  the  left  hand  side  is  positive.  However 

R 


we  can  expect  the  two  roots  away  from  xn+^  s  j  to  be  rather  large  since 
we  assume  3  is  small.  Thus  in  our  algorithm  to  find  the  root  we  select 
the  one  nearest  to  R/B. 

Note  that  for  nonlinearities  which  are  more  general  than 
cubic  a  Newton-Ralphson  scheme  could  be  used  on  each  step  with  first 
estimate  equal  to  R/B. 

3.3  A  Higher  Order  Method 

As  shown  above  Houbolt’s  implicit  method  is  accurate  to 


0(At  )  on  each  step  as  is  the  explicit  central  difference  method.  We 


now  seek  an  implicit  formula  with  errors  of  higher  order  than  AtH. 


This  can  be  done  by  simply  using  terms  in  a  backward  difference 
approximation  at  time  t  +  At;  the  error  will  decrease  as  we  use  more 
terms  in  the  difference  expression.  To  keep  computations  manageable 
we  investigated  up  to  ninth  order  schemes  and,  as  shown  later,  found 
the  eighth  order  scheme  the  most  stable.  This  scheme  replaces 
derivatives  by 


(13l32x  . -56l96x  +110754x  -132860x  ,+103320x  _-50652x  , 

n+1  n  n-1  n-2  n-3  n-4 


n+1 


2520  At' 


+I4266x  c-l764x  ,)  , 

- 5=5 - -HHit6) 


x  = 
n 


(13068x  . -35280x  +52920x  -58800x  „+44l00x  ,-21168x  . 

n+1  n  n-1  n-2  n-3  n-4 

5040  At 


+5880x  c-720x  ,  )  n 

- 2=5 - E±.  «<s 


On  substituting  into  (1)  we  obtain  . 


8x  +  B'x  .  =  R’ 
n+1  n+1 


which  can  be  solved  as  before  for  xn+^.  The  accuracy  of  the 

8, 


integration  is  0(At  )  on  each  step. 

A  comparison  of  Houbolt's  and  the  eighth  order  scheme  is 
presented  in  section  4.0 
3.4  Starting  Procedures 

Houbolt's  scheme  requires  x(t-2At),  x(t-At)  and  x(t)  in  order 
to  determine  x(t+At).  Thus  at  time  t  =  0  a  special  starting  procedure 
is  needed.  A  Taylor  series  scheme  is  an  obvious  choice.  Firstly,  on 
writing  the  equation  at  t  -  0  as 

(8) 

we  can  find  Xq  from  the  initial  conditions  Xq  =  p  and  xft  ■  q .  Then  x 


x0  +  cx0  +  kx0  +  exo  "  0 


0 


-1 


can  be  determined  from  the  Taylor  series 


and  x^  from 


x_i  =  xQ  -  At  xQ  +  -y  xQ  +  0(At  ) 


xx  =  xQ  +  AtxQ  +  ~~  Xq  +  0(At3) 


(9) 


(10) 


For  the  next  step  we  can  now  use  Houbolt's  scheme  since  we  know  x_^ ,  Xq 

4 

and  x^.  The  accuracy  of  Houbolt's  scheme  is  0(At  )  on  each  step. 

However  the  application  of  the  expansions  (9-10)  limits  the  accuracy  to 
3 

0(At  ).  Note  that  the  above  scheme  is  equivalent  to  some  authors' 
application  of  (9)  followed  by  a  central  difference  approximation  to 
(8)  about  t  =  0. 

g 

In  our  scheme,  which  is  accurate  to  0(At  )  on  each  step,  we 

require  a  knowledge  of  x(t-6At),  x(t-5At) .  x(t-At)  and  x(t)  in 

order  to  advance  to  x(t+At).  Thus  we  seek  a  starting  scheme  which 
determines  x_^,  x_2»  x_^»  xi »  x2  an<*  *3  to  an  accuracy  O(At^).  Our 
total  integration  scheme  will  then  be  accurate  to  O(At^).  Note  that  a 

starting  accuracy  higher  than  O(At^)  is  not  essential  as  in  one  cycle 

our  scheme  will  yield  an  error 

N  •  0(At8) 

2tf  7 

where  N  =  and  so  the  error  per  cycle  is  —  0(At  ).  We  employ  a 

Taylor  series  expansion  around  t  =  0  in  order  to  compute  x_^,  x_2>  x_^, 
x^,  X2  and  x^  i.e. 

.  (nAt)^  ”  ,  (nAt)3  —  ,  .  (nAt)8  VI  ,  rt,,^7x 

xn  "  x0  +  nAtx0  +  ~2!  x0  +  3~j  x0  +  •••  +  —5 ~  x0  +  0(At  > 


for  n  =  -3,  -2,  -1,  1,  2,  3,  where  Xq  is  found  from  substituting  the 

...  iv 

initial  conditions  into  the  differential  equation.  Xq,  Xq  and  higher 


derivatives  are  found  from  continued  differentiation,  viz 
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2» 


Kq  +  cXq  +  kxg  +  33xqXq  =  F  a)  cos  wt 


determines  Xq  while 

IV  ■'  2"  •  2  2 

xQ  +  cxQ  +  kxQ  +  30XqXq  +  S&XqXq  ■  -  F  o»  sin  wt 

IV 

determines  x^  ,  etc. 

Thus  we  have  a  starting  procedure  of  accuracy  O(At^)  which  is 
consistent  with  our  main  integration  scheme. 


3,5  Stability 

In  order  to  test  stability  of  schemes  of  higher  order  than 
Houbolt's  we  follow  [1]  to  derive  the  growth  matrix  on  each  step 
assuming  zero  damping  (c=0).  The  spectral  radius  (p)  of  the  matrix 
indicates  that  the  method  is  stable  if  p  <  1.  This  quantity  is 
computed  numerically  using  the  IMSL18  package  EIGRF  and  is  plotted  in 

9 

Fig.  2  for  Houbolt's  scheme  and  higher  order  schemes  up  to  0(At  ).  It 
can  be  seen  that  Houbolt's  scheme  is  unconditionally  stable  as  expected 
whereas,  for  example,  the  O(At^)  scheme  is  unconditionally  unstable  and 

g 

the  0(At  )  scheme  is  stable  provided  At/T  <  0.11  i.e.  more  than  10 
steps/cycle  are  required.  Since  the  latter  condition  Is  not  a 
practical  restriction  and  can  easily  be  imposed  we  will  use  this  scheme 
for  our  integration. 

The  higher  order  scheme  requires  more  computations  per  step 
than  Houbolt's  by  a  factor  of  about  2.  However,  since  a  much  larger 
step  length  is  possible  for  the  same  accuracy,  we  can  expect  an  overall 
increase  in  efficiency.  The  results  section  to  follow  shows  a 
comparison  of  efficiency  to  acquire  a  certain  accuracy.  Good  accuracy 
is  important  when  we  come  to  study  the  nonlinear  problem,  in  which  case 
the  starting  conditions  will  be  critical  in  deciding  which 
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amplitude-frequency  branch  the  solution  will  settle  upon.  Thus  the 
effect  of  the  starting  conditions  must  be  computed  accurately  as  we 
proceed  forward  in  time. 


4.0  RESULTS 

4.1  Linear  Equation 

Our  first  test  was  to  compare  our  eighth  order  results  with 
results  from  Houbolt's  scheme.  The  test  was  made  on  a  case  having 
higher  frequency  transients  i.e. 

x  +  O.lx  +  25x  *  25  sin  t 

with  Xq  =  0.1  and  x^  =  0.  It  can  be  seen  from  Figs.  3A  and  3B  that 
with  64  steps/cycle  the  eighth  order  scheme  gives  a  good  prediction  of 
the  exact  solution  whereas  Houbolt's  scheme,  as  expected,  has  large 
errors.  Even  with  256  steps/cycle,  Figs.  3C  and  3D,  Houbolt's  scheme 
is  showing  larger  errors  than  the  64  steps/cycle  eighth  order  scheme. 
Thus,  although  the  new  scheme  requires  twice  as  many  computations  per 
step,  the  number  of  steps  can  be  reduced  so  that  the  overall  efficiency 
is  about  twice  that  of  the  Houbolt  scheme  for  the  same  or  better 


accuracy. 


The  next  case,  Fig.  4,  is  the  simplest  case  of  no  damping 


x  +  x  ■  0 

It  can  be  seen  from  Fig.  4A  that  the  new  scheme  after  16  cycles  with  32 
steps/cycle  is  almost  exact  whilst  Houbolt's  scheme  shows  a  marked 
phase  shift  and  amplitude  decrement.  Even  with  128  steps/cycle  and 
after  16  cycles  (Fig.  4B),  a  slight  phase  shift  is  still  detected. 


Thus  the  new  scheme  is  shown  to  be  more  efficient  than 


Houbolt's  by  at  least  a  factor  2.  In  order  to  check  linear  accuracy  in 
other  cases  a  numerical  integration  of  both  the  linear  and  nonlinear 
equations  was  simultaneously  carried  out  in  all  our  further 
computations  and  the  linear  solution  was  always  compared  to  the  exact. 
Some  of  these  examples  are  covered  in  the  next  section. 

4.2  Nonlinear  Equation 

To  study  the  numerical  behaviour  in  the  nonlinear  case  we 
choose  for  our  model  the  constants 

c  =  0.1 

k  =  1  (11) 

0  *  0.1 

F  =  0.4 

The  response  curve  | A |  versus  u)  is  shown  in  Fig.  5A  for  0.2  <  u>  <  1.4 
and  in  Fig.  5B  for  1.18  <  w  <  1.38.  In  the  latter  case  we  observe  more 
clearly  bending  over  of  the  response  curve  to  the  right.  It  can  be 
seen  that  the  Kryloff  and  Bogoliuboff  approximation  given  earlier 
(Eq.  5)  is  much  more  accurate  than  the  more  frequently  used  formula 
(Eq.  3)  in  the  peak  region  A  *  3.  Also  shown  is  Duffing's 
approximation  based  on  zero  damping  (c»*0). 

Our  experience  in  generating  the  response  curve  showed  that 
steady  state  solutions  were  obtained  in  typically  10-40  cycles  of 
period  2tt/w.  Away  from  the  nonunique  parts  of  the  curve  i.e.  u>  <  1.18 
and  w  >  1.31  solutions  were  obtained  in  fewer  iterations.  In  the  range 

1.18  <  u  <  1.31  the  branch  of  the  response  on  which  the  solution  lay 

•  • 

depended  upon  the  initial  estimates  xn,  xn.  For  instance  with  xn  *  0 
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several  computations  were  made  with  various  starting  values  x^.  With 
Xq  =  2  the  jump  from  the  lower  branch  occurred  at  about  w  =  1.22  while 
with  Xq  =  -3.5  the  jump  occurred  at  about  oj  =  1.32  as  shown  on  Fig. 

5B. 

Other  starting  values  produced  jumps  as  shown  on  Fig.  5C, 
which  also  shows  a  finer  determination  of  the  cut  off  jumps  for  non 
uniqueness.  The  corresponding  phase  angles  relative  to  the  driving 
force  are  given  in  Table  1.  It  was  initially  expected  that  there  would 
be  a  trend  of  the  jump  values  of  ui  with  changing  Xq.  However  as  can  be 
seen  there  is  no  apparent  trend  in  the  jump-rather  there  seems  to  be  a 
randomness  which  is  unexplainable.  This  is  disconcerting  particularly 
as  there  is  no  consistent  way  to  produce  the  limits  of  the  upper  branch 
to  the  right  and  of  the  lower  branch  to  the  left.  These  limits  in  our 
case  were  determined  by  trial  and  error  to  be  near  w  *  1.18  on  the 
lower  branch  and  near  =  1.31  on  the  upper  branch  taking  x^  =  4 
for  the  lower  branch  and  Xq  =  -3.5  for  the  upper  branch.  In  these 
areas  near  the  vertical  slopes  solutions  would  often  take  a  long  time 
to  settle  on  the  final  branch.  The  amplitude  would  tend  to  stick 
around  the  value  on  one  branch  for  many  cycles  and  then  quite  suddenly 
the  amplitude  would  change  to  a  value  on  the  other  branch.  An  example 
of  this  is  shown  in  Fig.  6  for  Xq  *  -3.5  and  u>  “  1.32. 

The  unstable  part  of  the  response  curve  shown  as  path  3-6  on 
Fig.  1  was  never  predicted  numerically  and  did  not  cause  any  problems. 
The  final  amplitude  would  either  lie  on  the  lower  or  upper  stable 


branches 
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The  drastic  change  in  solution  from  one  branch  to  the  other 
can  be  seen  by  comparing  Fig.  7A  to  Fig.  7B.  These  figures  show 
converged  solutions  over  one  period  of  the  driving  force  i.e.  2ti / to  for 
u)  =  1.20  and  oo  =  1.22  respectively  (with  Xq  =  2).  It  can  be  seen  that 
there  is  a  large  phase  shift  as  well  as  an  amplitude  change. 

Since  the  branch  solution  chosen  by  the  numerical  procedure 
is  so  dependent  on  initial  conditions  it  was  considered  worth 
investigating  the  effect  of  step  size  At  on  the  branch  chosen.  This  is 
because  a  numerical  truncation  error  might  affect  the  solution  to  such 
an  extent  that  it  will  jump  at  a  different  u>  value  for  the  same  x^. 

Thus  we  repeated  the  computations  of  Fig.  5C,  which  use  64  steps/cycle, 
with  128  steps/cycle.  It  was  found  that  the  results  were  identical, 
except  for  x^  =  10,  when  the  jump  shifted  right  over  to  between  1.30 
and  1.32  (compared  to  1.18  to  1.20).  However  there  was  still  no  trend 
in  results  even  with  the  smaller  truncation  error. 

5.0  CONCLUSIONS 

As  a  first  step  to  studying  the  nonlinear  response  of 
structural  components  we  have  investigated  numerically  a  one-degree-of- 
freedom  nonlinear  response  equation.  It  has  been  shown  that  the 
multivalued  response  curve  of  amplitude  versus  frequency  can  be 
obtained  using  a  numerical  method  of  integration  which  is  more 
efficient  than  the  more  standard  Houbolt's  scheme.  The  jump  phenomenon 
from  one  branch  of  the  response  curve  to  another  has  been  investigated 
and  the  branch  chosen  is  found  to  depend  upon  the  initial  conditions 
hut  no  trend  can  be  detected  in  relating  the  frequency  at  which  the 
jump  occurs  with  initial  conditions. 
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Table  1 .  Phase  Angle,  $,  in  degrees. 
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FIG.  1:  SCHEMATIC  OF  THE  RESPONSE  CURVE  FOR  0  >  0.  SHOWING 

JUMP  PHENOMENA 


FIG.  2:  STABILITY  OF  VARIOUS  DIFFERENCE  SCHEMES  FOR  X 
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FIG.  3(a):  ACCURACY  OF  THE  EIGHTH  ORDER  SCHEME  WITH  64  STEPS/CYCLE 
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FIG.  3(b):  ACCURACY  OF  THE  EIGHTH  ORDER  SCHEME  WITH  64  STEPS/CYCLE 
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FIG.  3(d):  ACCURACY  OF  HOUBOLT'S  SCHEME  WITH  256  STEPS/CYCLE 
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FIG.  4(a):  ACCURACY  OF  EIGHT  ORDER  SCHEME  FOR  SIMPLE  HARMONIC 

WITH  32  STEPS/CYCLE 
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FIG.  4(b):  ACCURACY  OF  HOUBOLT'S  SCHEME  FOR  SIMPLE  HARMONIC  MOTION 

WITH  128  STEPS/CYCLE 
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FIG.  7(a):  SUDDEN  CHANGE  IN  SOLUTION  ACROSS  THE  JUMP 
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FIG.  7(b):  SUDDEN  CHANGE  IN  SOLUTION  ACROSS  THE  JUMP 
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